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ABSTRACT 

A  barotropic  primitive  equation  model  on  a  global  tropi- 
cal grid  using  operational  real  data  is  used  to  test  various 
boundary  conditions  and  methods  of  initialization  including 
numerical  variational  analysis.   Comparisons  of  forecast 
accuracy  are  made  between  staggered  and  non-staggered  grids. 

All  prediction  models  produce  better  verification  stat- 
istics than  persistence,  with  the  staggered  grid  verifying 
better  than  the  full  grid.   This  appears  to  be  a  result  of 
two-gridlength  noise  being  of  greater  magnitude  in  the  full 
grid  model.   Less  altering  of  individual  synoptic  systems 
occurs  in  the  variational  analysis  initialization  compared 
to  two  forms  of  the  nonlinear  balance  equation.   As  a  re- 
sult, nondivergent  variational  analysis  input  fields  to  the 
prediction  models  result  in  the  best  verification  values. 
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I.   INTRODUCTION 

Numerical  weather  prediction  in  the  tropics  has  devel- 
oped slower  than  prediction  in  the  mid-latitudes  for  many 
reasons.   The  lack  of  data  and  less  understanding  of  the 
actual  physical  processes  governing  the  circulations  in  the 
tropics  are  probably  the  major  deterrents  to  progress.   In 
recent  years,  the  increase  of  meteorological  data  and  in- 
creased understanding  of  the  dynamics  in  the  tropics  have 
sparked  interest  in  attempts  at  operational  forecasting. 

Tropical  weather  prediction  can  be  approached  on  either 
a  global  band  grid  or  a  global  grid.   Certain  characteris- 
tics of  tropical  prediction  make  the  global  band  advanta- 
geous over  the  global  grid.   In  mid-latitudes,  much  of  the 
baroclinic  energy  exchange  occurs  on  the  long  wave  scale 
while  in  the  tropics,  the  input  of  energy  on  the  cumulus 
scale  is  also  an  important  contributor  to  the  dynamics. 
This  complicates  the  global  prediction  by  making  the  tro- 
pics a  subset  of  the  entire  globe.   In  mid-latitude  regions 
of  the  Northern  Hemisphere,  data  are  routinely  available  at 
all  the  mandatory  levels  as  well  as  significant  levels. 
In  the  tropics,  data  availability  is  a  different  story. 
Upper  level  data  near  300  mb  are  available  from  aircraft 
reports  combined  with  scarce  conventional  rawinsondes. 
Gradient  level  and  surface  level  data  are  also  routinely 
available  in  numbers  approaching  those  necessary  for  input 


to  an  operational  model.   However,  the  tropical  regions  are 
lacking  in  mid- tropospheric  observations.   The  combination 
of  a  smaller  scale  of  motion  and  the  observational  data 
available  make  data  initialization  in  the  tropics  somewhat 
different  than  that  in  mid-latitudes.   Reduction  of  the 
divergence  equation  to  the  balance  equation,  which  is  a 
generalized  diagnostic  equation  used  to  relate  wind  and 
heights  in  mid-latitudes,  is  dependent  upon  a  scale  analysis 
in  which  the  Rossby  number  is  assumed  less  than  one.   This 
is  not  satisfied  in  the  tropics  however.   Also,  the  optimum 
amount  of  divergence  remaining  in  the  initialized  fields, 
which  are  used  as  input  to  a  prediction  model,  may  differ 
in  the  tropics  from  that  in  mid-latitudes. 

As  a  result  of  the  above  considerations,  Fleet  Numeri- 
cal Weather  Central,  Monterey  now  produces  an  operational 
wind  and  temperature  analysis  for  a  global  band  from  41S  to 
60N  at  eight  levels  ranging  from  the  surface  to  100  mb .   A 
more  detailed  description  is  presented  later  in  this  thesis. 
This  analysis  is  a  potential  input  to  the  primitive  equation 
prediction  model.   Steinbruck  (1971)  conducted  preliminary 
experiments  on  this  global  band  grid  using  climatological 
data.   He  investigated  the  properties  of  various  boundary 
conditions  on  the  primitive  equation  prediction  model  and 
various  diagnositic  balancing  schemes,  including  numerical 
variational  analysis  (NVA)  developed  by  Dr.  J.  M.  Lewis 
(1972)  for  the  global  band. 

The  purpose  of  this  thesis  is  to  extend  this  research 


to  daily  operational  real  data  and  to  introduce  a  staggered 
grid  system  in  an  attempt  to  reduce  memory  and  time  require- 
ments on  the  computer.   A  deliberate  attempt  is  made  to 
avoid  introducing  averaging  and  diffusion  in  order  to  pro- 
duce as  pure  a  forecast  as  possible  and  isolate  true  differ- 
ences between  boundary  conditions  and  the  various  techniques 
for  calculating  initial  data. 
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II.   PRIMITIVE  EQUATION  PREDICTION  MODEL 


A.   BASIC  MODEL 

A  ten-layer  primitive  equation  model,  developed  by  Har- 
rison (1970)  and  Elsberry  and  Harrison  (1970,  1972),  was 
used  as  a  basis  for  the  barotropic  models  used  in  these  ex- 
periments.  The  equations  for  the  ten-layer  model  are: 


i£  =  -L(u)  +  fv  -  mf£  +  Fx 

d  t  dx 

lX  =  -L(v)  -  fu  -  mil  +  Fy 
3t  3y 

||  =  -L(6)  +  heat 


ia   =  - 


3t 
3<J> 


-L(q)  +  moisture 


1000 
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P     v1000' 
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T  t  ~\  2r3    us     3    vs 

L(s)  =  m   [- —  / — \  +  - —  /_\i 
l3x  ^m  >         3y  U  'J 


iy  vm  'J  '  3p  (ws) 
and  normal  meteorological  symbols  are  used  with 
A  =  finite  difference 
=  layer  mean  value 
1000  =  1000  mb  surface 
m  =  map  factor 
S  =  general  scalar  dependent  variable 


(1) 
(2) 

(3) 

(A) 

(5) 
(6) 

(7) 
(8) 
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L  =  advective  operator 
Linear  stability  is  maintained  by  requiring  the  ratio 
of  the  space  increment  to  the  time  increment  be  greater 
than  the  phase  speed  of  the  fastest  waves,  which  are  gravity 
waves.   At  300  mb ,  and  on  the  grid  used,  this  time  step  is 
five  minutes  since  the  space  increment  near  60N  is  approxi- 
mately 75  nautical  miles. 

All  models  use  a  forward  time  step  in  the  first  predic- 

1     0 
tion  step:   S   =  S   +  A  t  *  function  (S°).   All  subsequent 

prediction  steps  use  a  leapfrog  scheme  where:  Sn    =  Sn 

+  2  At  *  function  (Sn)  . 

B.   NON-STAGGERED  GRID  BAROTROPIC  PREDICTION  MODEL 

Because  the  prediction  models  to  be  discussed  here  are 
barotropic,  Eqns  .   (3),  (4)  and  (7)  are  not  needed  and  since 
all  prediction  is  done  at  the  300  mb  level,  Eqn.  (5)  be- 
comes : 

8({)300       2   3    »<J>     3_  ,v± 
3t     =  -  m   [g^  (m  )  +  3y  (m  )] 

The  frictional  terms  in  Eqns.  (1)  and  (2)  are  neglected. 

Input  data  to  the  models  are  the  winds  and  geopotential 
at  300  mb . 

The  finite  difference  formulas  are  based  on  Arakawa's 
(1966)  flux  form.   This  form  of  differencing  attempts  to 
preclude  nonlinear  computational  instability  by  requiring 
that  the  advective  terms  conserve  the  mean  square  of  any 
advected  quantity.   Therefore,  mean  square  kinetic  energy 
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is  conserved  by  the  nonlinear  terms  and  nonlinear  instabili- 
ty is  avoided.   All  horizontal  derivatives  are  centered  two 
grid  mesh  derivatives.   The  horizontal  advective  terms  for 
a  general  variable  S  with  J  the  index  in  the  y  direction  and 
I  the  index  in  the  x  direction  are  written: 

8    Su  Su  Su 

9x  (m  )   _   =  <m  )  T  TJ,  ..   -  (  m) 

J, I         J, 1+1/2         J, 1-1/2 

Ax 


=  [  SJ,I+1+  SJ»I    1   u      +  u 

2        J  2  S,I+1  »J,I] 

Ax 
J»I  +  SJ,I-1   1   u      +  u 

[       2 3  2  [mJ,I     mJ,I-l] 

Ax 

Note  that  this  finite  differencing  only  averages  parameters 
over  two  grid  points  and  one  grid  length. 


C.   STAGGERED  GRID  BAROTROPIC  PREDICTION  MODEL 

Computation  of  space  derivatives  in  finite  difference 
form  incorporates  data  from  adjacent  grid  points.   This 
prompted  attempts  at  distributing  the  variables  at  alternat- 
ing grid  points  in  an  attempt  to  reduce  computer  time  and 
memory.   This  approach  was  first  attempted  by  Eliassen 
(1956)  and  subsequently  has  been  used  by  many  investigators 
in  various  forms.   The  prediction  and  conservation  equa- 
tions on  the  staggered  grid  are  the  same  as  on  the  full 
grid.   The  staggered  grid  distribution  used  in  the  baro- 
tropic  prediction  model  involves  u  and  v  winds  at  alternat- 
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ing  points,  (j)  and  to  at  diagonally  alternating  points  as 
shown  in  Fig.  1.   The  Mintz-Arakawa  general  circulation 
model  as  described  by  Langlois  and  Kwok  (1969)  and  an  ocean 
circulation  model  by  Haney  (1971)  also  use  this  form  of  the 
staggered  grid.   This  grid  uses  the  least  number  of  points 
possible  and  still  retains  the  characteristic  that  the 
pressure  gradient  force  is  evaluated  over  the  same  space 
increment  as  in  the  non-staggered  grid.   Advective  terms, 
however,  involve  a  larger  space  increment. 

The  horizontal  advective  terms  for  a  wind  variable  X 
are  written: 

|_  [XU]   =  (XU)   _    (XU) 
dx   m        m  m 

x=J,I    x=J, 1+1/2   x=J, 1-1/2 

2Ax 


[Xj'I+1+  ^3i  [(S)         +(-)      ]  - 

2         J2  L  Vj.I+1   m  J, I 


2Ax 


XT  T  +  X 


[  J>1  ,  J>1"13  \   [(  *)    +  ft    ] 
I 2    m  j, i    mj,i-i 

2Ax 

Refer  to  Fig.  1  for  indexing.   In  a  similar  manner,  the 
horizontal  advection  of  geopotential  (4>)  takes  the  form 

3x   <-m  >       ~    ^m  >  Cm  > 

0=J,I    q=j, 1+1/2     0=J, 1-1/2 

2Ax 

..^,1+1  + V   i    *       +(|         , 

2 2    m  J, I     m  Jt1,I 

2Ax 
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FORECASTING  GRIDS 
GLOBAL   BAND:      49   x   146    (COLUMN   145   -   1,    COLUMN  146   =■   2) 
STAGGERED  GRID:         U,V      -   25   x   74    (COLUMN  73   =    1,    COLUMN  74 
i.aj     -  26   x   74    (COLUMN  73   =   1,    COLUMN   74 
STAGGERED  DATA  EXTRACTION 
U.V   (J-1,25,      1=2,73)      =     U.V    (2*J-1,    2*1-1) 
*.u>   U-2,25,      1=2,73)      =     J.O)   (2*J-2,   2*1-2) 
ROW  1   AND  26   DETERMINED   BY  BOUNDARY  CONDITIONS. 


Fig.  1.   Staggered  grid  description 
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_  r*J»I+  ^,1-1,    1  r     ,X,  ,  (X, 

1  2  3    2  [  gj>M+  gj=i!i-i3 

2Ax 

The  vertical  motion  used  in  vertical  advection  of  wind 

speed  is  obtained  by  averaging  the  four  nearest  diagonal 

values  of  a) . 

As  in  the  full  grid,  this  form  of  finite  differencing 
averages  parameters  over  only  two  points  but  results  in 
differencing  over  two  full  grid  lengths. 

Vertical  motion  is  computed  from  the  continuity  equa- 
tion using  winds  that  are  obtained  by  averaging  predicted 
winds  over  two  grid  increments  located  diagonally  from  the 
vertical  motion  gridpoint. 

D.   BOUNDARY  CONDITIONS 

The  computationally  stable  integration  of  the  primitive 
equations  requires  a  set  of  boundary  conditions  that  are 
consistent  with  the  finite  difference  scheme. 

Flux  through  the  top  of  the  domain  is  eliminated  by  set- 
ting w  =  0  at  100  mb . 

Two  sets  of  horizontal  boundary  conditions  were  tested 
in  the  prediction  model: 

1 .   No-flux  boundaries 

Elsberry  and  Harrison  (1971)  discuss  the  boundary 
conditions  that  were  applicable  to  the  basic  ten-level  model 
and  which  are  used  in  one  version  of  the  barotropic  model. 
In  the  full  grid,  no  flux  of  mass  or  energy  is  allowed  at 
the  northern  boundary  by  setting: 
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m 


J+l 


J+l 


m 


VT  ,  u 


J+l 


=  u 


^J+l=  *. 


where  J+l  is  the  northern-most  row  in  the  grid.   The  south- 
ern boundary  is  treated  in  a  similar  manner.   These  boundary 
conditions  result  in  the  model  boundary  being  a  fictitious 
grid  row  between  the  two  outer  rows  of  input  data  across 
which  there  is  no  flux  of  mass  or  energy.   East-west  bound- 
ary conditions  are  not  needed  since  the  global  band  is  con- 
tinuous . 

The  u  boundary  conditions  on  the  staggered  grid  are 
similar  to  those  on  the  full  grid.   In  order  to  conserve 

kinetic  energy  on  the  staggered  grid  <J)     =  |   where  J+l  is 

J+ 1     J 

the  northern-most  row  of  the  staggered  geopotential  values. 
To  satisfy  the  no-flux  criterion,  the  outer  row  of  v  com- 
ponent wind  must  be  identically  zero.   This  specification 
of  boundary  conditions  results  in  the  model  boundary  being 
coincident  with  the  outer-most  row  of  staggered  grid  wind 
components . 

2 .   Restoration  boundaries 

Kesel  and  Winninghoff  (1970)  describe  a  constant 
flux,  restoration  boundary  technique.   This  technique  is 
used  with  initial  data  from  the  nonlinear  balance  equation 
which  has  consistent  boundary  conditions,  and  balancing 
produced  by  a  numerical  variational  analysis  scheme.   Both 
of  these  techniques  are  later  described  in  detail.   The 
restoration  technique  involves  restoring  the  newly  predicted 
values  near  the  boundaries  toward  their  value  at  the  previous 
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time  step  with  a  specified  restoration  coefficient.   In  the 
full  grid  prediction,  these  coefficients  are  a  linear  varia- 
tion of  1.0,  0.8,  0.6,  0.4,  0.2,  0.0.   This  results  in  a 
persistence  forecast  on  the  outer  boundary,  pure  dynamical 
forecast  six  rows  in  and  a  linear  combination  between.   Note 
that  the  data  boundary  is  the  model  boundary.   In  the 
staggered  grid  a  similar  approach  is  used  with  coefficients 
of  1.0,  0.6667,  0.333,  0.0.   In  the  staggered  grid,  a 
fictitious  geopotential  value  is  needed  one  grid  length 
outside  the  boundary  of  the  available  data  so  a  linear 
extrapolation  on  the  input  data  geopotential  gradient  is 
used  : 


<!> 


J+2,I 


=    2    <t>J+1    -    cf> 


This  results  in  the  model  geopotential  boundary  being  one 
row  outside  the  input  data  boundary  and  the  model  wind 
boundary  coinciding  with  the  data  boundary. 
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III.   DATA  DESCRIPTION 

First  test  experiments  with  the  prognostic  and  diagnos- 
tic models  used  a  climatological  input  produced  by  Jenne 
(1969)  of  the  National  Center  for  Atmospheric  Research.   In 
particular,  these  experiments  used  300  mb  January  climato- 
logy for  the  area  of  the  global  band.   Although  this  data 
lacks  the  small-scale  features  which  are  important  in  day- 
to-day  forecasting,  the  strong  interaction  between  the  tro- 
pics   and  the  middle  latitudes  through  the  subtropical  jet 
is  present.   In  addition,  boundary  condition  effects  and 
development  of  near-equatorial  circulations  can  be  tested. 
The  300  mb  level  was  chosen  so  that  comparable  later  experi- 
ments using  real  data  would  allow  a  maximum  of  actual 
observations . 

The  real  data  utilized  in  later  experiments  are  from 
the  Fleet  Numerical  Weather  Central  global  band  analysis 
developed  by  Grayson  (1971).   This  analysis  covers  the 
entire  global  band  from  60N  to  41S.   The  grid  dimensions 
are  49  X  144  on  a  Mercator  secant  projection  map  true  at 
22.5°  latitude.   The  mesh  length  is  2.5°  longitude  (=150 
n.mi.)  at  the  equator  and  reducing  to  75  n.mi.  at  60N.  This 
grid  size  has  the  capability  of  depicting  relatively  small 
scale  synoptic  events,  provided  the  observations  are  suffi- 
ciently dense.   Grayson  (1971)  used  a  modified  version  of 
the  successive  approximations  technique  to  analyze  sea 
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level  pressure  and  surface  wind.   This  approach  has  subse- 
quently been  extended  to  the  analysis  of  seven  levels  of 
upper  level  winds,  but  it  should  be  noted  that  no  vertical 
consistency  check  is  applied.   In  an  attempt  to  accommodate 
the  large  variation  of  data  coverage  and  the  various  scales 
of  motion  in  the  wind  analysis,  the  area  influenced  by  each 
observation,  in  addition  to  being  dependent  on  the  scan 
number,  is  dependent  upon  the  reported  wind  direction  (if 
greater  than  13  knots).   This  approach  was  suggested  by 
Endlich  and  Mancuso  (1968)  . 

The  previous  six  hour  old  analysis  is  used  as  a  first 
guess  field.   450  reports  per  analysis  is  a  typical  number 
of  reports  at  300  mb  with  approximately  the  following  lati- 
tudinal distribution: 

45N  -  60N  =  37% 

30N  -  45N  =  33% 

15N  -  30N  =  20% 

15S  -  15N  =   6% 

30S  -  15S  =   2% 

41S  -  30S  =   2% 
No  upper  air  geopotential  fields  are  analysed  by  the 
global  band  analysis  program.   Therefore,  a  first  guess 
field  is  produced  using  a  combination  of  the  current  FNWC 
hemispheric  analysis  north  of  27N  and  climatology  south  of 
17N  with  a  linear  combination  in  between. 

Since  the  real  winds  and  climatological  heights  along 
the  southern  boundary  of  the  global  grid  are  not  geostrophi' 
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cally  consistent,  the  height  values  are  altered  along  the 
boundary.   Geostrophic  gradients  appropriate  to  the 
analysed  v  component  of  the  wind  are  added  to  the  average 
climatological  height  along  the  southern  wall.   This 
approach  is  applied  to  all  diagnostic  models  except  no-flux 
boundary  conditions  where  it  is  not  applicable. 
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IV.   INITIALIZATION  METHODS 

A.   NONLINEAR  BALANCE  EQUATION  WITH  NO-FLUX  BOUNDARY 
CONDITIONS 

Elsberry  and  Harrison  (1971)  utilized  the  nonlinear 
balance  equation  and  applied  no-flux  boundary  conditions  to 
formulate  a  diagnostic  model  which  matches  the  no-flux  ver- 
sion of  the  prediction  model.   The  nondivergent  part  of  the 

2 
wind  is  obtained  from  £  =  V  \p   where  £  is  the  relative  vorti- 

city  obtained  from  the  observed  u  and  v  wind  components. 

No-flux  boundary  conditions  imply  that  the  north  and  south 

boundaries  are  streamlines  with  no  v  component  wind  across 

them.   Conservation  of  kinetic  energy  requires  that  the 

two  outer  rows  of  <J)  at  the  boundaries  be  set  equal.   In 

order  to  avoid  setting  up  gravity  waves  propagating  along 

the  boundary,  the  geopotential  along  the  entire  boundary 

is  set  equal  to  the  average  value.  \b         .,  is  set  equal  to 
n  or  soutn  ^ 

ip    .,  -uAy  where  u  is  the  mean  zonal  wind  component 
Ynorth     J  r 

averaged  over  the  grid  and  Ay  is  the  distance  from  the  north 
to  south  boundary.   The  natural  cyclic  east-west  boundary 
conditions  alleviate  the  need  for  artificial  boundaries  in 
this  direction . 

The  nondivergent  wind  is  calculated  from  the  relaxed  ty 
field  using  v  =  k  x  V^  .   No-flux  conditions  are  forced  by 
making  v(49,I)  =  v(48,I)  m(49)/m(48)  and  similarly  at  the 
southern  wall  where  m  is  the  map  factor.   Then,  the  geo- 
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,  2 

potential  field  is  generated  by  relaxing  V  (J)  using  the  non- 
linear balance  equation: 

V    <j>    =    Vf    •    V>    +    fV    \\>    +    2J    (f-&     |^-) 

dx       dy' 

One  should  note  that  large  deviations  in  the  geopotential 
from  the  input  geopotential  are  to  be  expected  due  to  the 
zonal  averaging  of  the  boundary  values. 

B.   NONLINEAR  BALANCE  EQUATION  WITH  RESTORATION  BOUNDARY 
CONDITIONS 

Kesel  and  Winninghoff  (1970)  describe  a  restoration 
technique  which  eliminates  the  need  for  altering  the  bound- 
ary geopotential  field  and  appears  to  control  false  reflec- 
tion of  physical  and  computational  modes  at  the  boundary. 
Wind  and  height  values  of  the  input  data  are  not  altered 
on  the  north-south  boundaries  in  contrast  to  the  no-flux 
method.   The  technique  involves  restoring  after  each  time 
step  the  values  of  the  parameters  near  the  walls  toward 
the  value  in  the  previous  time  step  with  a  specified  restor- 
ation coefficient.   This  results  in  a  purely  dynamic  fore- 
cast in  the  interior,  a  persistence  forecast  on  the  boundary 
and  a  combination  in-between.   The  technique  results  in  the 
boundaries  acting  as  an  energy  sponge   for  externally  pro- 
pagating meteorological  and  gravity  waves. 

When  working  with  climatological  fields  for  input  data, 
the  southern  boundary  of  ip  was  determined  similar  to  the  no- 
flux  case  by  using  a  column  average  of  the  zonal  wind  com- 
ponent.  This  technique  applied  to  real  data,  however,  pro- 
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duced  unrealistic  meridional  velocity  components  along  the 
southern  wall  due  to  the  large  and  rapid  variations  in  the 
east-west  direction  of  the  average  zonal  wind  component. 
To  correct  this  problem,  the  ip's  calculated  as  above  are 
averaged.   The  southern  ty    value  is  set  to  this  value  which 
results  in  a  north-south  \p    gradient  consistent  with  the 
mean  zonal  wind.   Upon  this  average  value  at  the  southern 
boundary,  variations  are  superimposed  to  reflect  east-west 
gradients  which  are  consistent  with  the  analysed  v  compon- 
ents on  the  southern  wall.   The  remainder  of  the  balancing 
is  identical  to  the  no-flux  approach. 

C.   NUMERICAL  VARIATIONAL  ANALYSIS 

Lewis  (1972)  has  adapted  Sasaki's  (1958)  variational 
method  to  obtain  dynamically  consistent  initial  fields  from 
input  objective  operational  analyses.   The  fields  are  ad- 
justed so  they  satisfy  the  equations  of  motion  and  minimize 
the  functional: 
I  =  /  /   [S(u-u)   +  a(v-v)   +  B(Z-Z)   +  a(|^-)+  a(|^-)]8s  (9) 

with  u,  v,  z  being  the  input  analysis.   a,  3  are  weighting 
functions  specified  beforehand  which  can  be  varied  in  res- 
ponse to  the  density  and  reliability  of  the  input  observa- 
tions and  which  determine  how  closely  the  desired  fields  u, 
v,  and  z  fit  the  input  analysis.   The  dynamic  weight  a  is  a 
similar  function  which  controls  the  horizontal  acceleration 
(calculated  from  the  momentum  equations  given  below)  and 
directly  couples  the  wind  and  height  information  via  the 
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assumed  dynamical  relations. 

3u       3Z    —    T7   _i_  r  3v        8Z    —    T7     c 

-|—  =  -  g- v  •  Vu  +  fv  —_  «  -   g^-  "  v  •  Vv  -  fu 

dt       8x  dt       &  dy 

The  use  of  these  equations  as  dynamical  constraints  require 

the  simultaneous  solution  of  the  three  equations  for  u,  v, 

and  z  (so-called  Euler  equations)  .   To  simplify  the  Euler 

equations,  the  inertial  terms  in  the  dynamical  constraints 

are  modified : 

lii  =  -  eM  -  v  •  Vu  +  fv        IZ  m    -  g3jZ  _  v  .  y^  _  fu 
3t      88x  dt       3y 

As  in  the  restoration  balancing,  the  climatological  southern 

height  values  are  averaged  and  altered  to  make  them  geo- 

strophically  consistent  with  the  v  component  winds  before 

the  NVA  analysis  is  made  on  the  fields. 
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V.   DIAGNOSTIC  RESULTS 

Diagnostic  experiments  were  aimed  at  three  objectives: 
compare  NVA  results  with  various  forms  of  the  nonlinear 
balance  equation;  determine  NVA  weighting  constants  which 
most  advantageously  initialize  the  prediction  model;  and 
test  various  methods  of  limiting  the  divergence  in  the  NVA 
output . 

The  first  objective  was  to  compare  the  results  from  the 
variational  analysis  approach  against  those  using  the  non- 
linear balance  equation  with  no-flux  and  restoration  bound- 
ary conditions.   One  should  first  note  the  basic  difference 
in  solution  approach  between  the  two  methods.   The  non- 
linear balance  alters  the  geopotential  field  based  on  the 
nondivergent  wind  input  from  the  streamf unction  solution, 
while  the  variational  analysis  simultaneously  alters  both 
the  winds  and  geopotential. 

Area  average  wind  changes  listed  in  the  tables,  and 
wind  changes  analysed  in  the  charts  to  follow  are  obtained 
by  subtracting  the  input  field  from  the  balanced  field. 
Therefore,  a  positive  Value  means  an  increase  of  westerly 
winds  In  the  u  component  and  an  increase  in  southerly  winds 
in  the  v  component.   Area  average  modulus  change  is  the 
average  of  the  absolute  value  of  the  individual  wind 
changes  which  gives  a  measure  of  the  magnitude  of  the 
average  change  at  each  grid  point.   The  values  given  are  an 
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average  of  the  u  and  v  component  changes  which  were  found 
to  have  similar  values.   The  RMS  change  is  the  square  root 
of  the  average  of  the  squares  of  the  balancing  changes. 
RMS  u  and  v  components  are  averaged  since  similar  values 
were  obtained. 

Fig.  2  shows  a  typical  v  wind  component  input  to  the 
objective  analysis  programs.   Fig.  3  shows  the  v  component 
changes  produced  by  the  nonlinear  balance  solution  with 
restoration  boundary  conditions.   The  nonlinear  balance 
approach  tends  to  reduce  the  absolute  magnitude  of  the  wind 
maxima  in  both  the  u  and  v  component  wind  with  the  centers 
of  maximum  change  being  of  nearly  equal  magnitude.   This 
wind  change  represents  the  divergent  component  of  the  wind. 
As  shown  in  Table  1,  the  modulus  and  root-mean-square  change 
in  the  no-flux  approach  are  larger  than  in  the  restoration 
approach  due  to  the  crude  boundary  approximations.   The 
non-realistic  boundary  conditions  in  the  no-flux  approach 
only  produce  significant  differences  from  the  restoration 
boundaries  within  about  five  rows  of  the  boundary.   The  wind 
speed  changes  produced  by  the  variational  analysis  approach, 
as  displayed  in  Fig.  4,  show  smaller  magnitudes  in  a  larger 
scale  pattern.   Generally,  the  tendency  is  to  reduce  the 
wind  maxima  and  increase  the  intervening  wind  minima.   Table 
1  shows  an  average  modulus  wind  change  of  0.60  m/sec  in  the 
pure  NVA  analysis  as  compared  to  a  modulus  change  approxi- 
mately four  times  as  large  for  the  nonlinear  balance  method. 
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Fig.  2.   300  mb  Global  band  v  component  wind  analysis  from 
60E  to  150W  for  08/1200Z  January  1972.   Units  = 
i/sec . 
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Fig.  3.   V  component  wind  speed  changes  resulting  from  non- 
linear balance  solution  with  restoration  boundary 
conditions.   Positive  (+)  value  indicates  increase 
in  southerly  wind  speed.   Units  =  m/sec. 
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TABLE  1 
METHODS  OF  BALANCING 


Balancing  Area     Area     RMS      Area     Area      RMS 
Type     Average  Average  Wind     Average  Average   Height 
Wind     Modulus  Change   Height   Modulus   Change 
Change   Wind     (m/sec)  Change   Height    (m) 
(i/sec)  Change           (m)      Change 
(m/sec) (m) 

I.  Resulting  Field  Nondivergent 

NVA 

(1,7,25)   u=.50    1.81     2.61     -22.0    31.0      37.0 

Made       v=.27 

Nondivergent 

No-flux    u=-.26   2.84     4.12      96.0   131.0     175.2 
Balance    v=  .23 
Equat  ion 

Res  toration 

Balance    u=-.03   1.89     2.87     106.7   131.0     174.4 

Equation   v=  .24 

II.  Resulting  Field  Divergent 

NVA        u=  .52    .60     1.10     -22.0    31.0      37.0 
(1,7,25)   v=  .025 
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Fig.  4.   V  component  wind  speed  changes  resulting  from  NVA 
initialization  program  (1,7,25).   Positive  (+) 
value  indicates  increase  in  southerly  wind  speed. 
Units  =  m/sec . 
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Fig.  5.   300  mb  Global  band  height  field  from  60E  to  150W 
for  08/1200Z  January  1972.   Units   =  decameters. 
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Fig.  5  shows  the  08/1200Z  January  1972  global  band 
height  field  used  in  these  experiments.   In  evaluating 
height  changes  produced  by  the  various  methods,  one  must 
remember  that  the  geopotential  field  was  obtained  from  the 
derived  nondivergent  wind  analysis.   In  addition,  approxi- 
mately two-thirds  of  the  height  field,  that  south  of  17N, 
is  a  climatology  field.   Keeping  this  in  mind,  one  can 
see  from  Table  1  that  the  NVA  analysis  alters  the  heights 
only  about  one-quarter  as  much  as  the  nonlinear  balance 
approach.   In  the  region  with  real  height  data,  the  NVA 
analysis  tends  to  fill  troughs  and  reduce  ridge  heights  in 
order  to  achieve  balance.   These  changes  appear  to  be  dis- 
placed toward  the  downwind  side  of  the  synoptic  feature  in 
the  mid-latitudes.   The  scarcity  of  data  in  the  southern 
hemisphere  tends  to  produce  isolated  wind  maxima  and  the 
NVA  introduces  height  gradients  to  partially  balance  these 
winds.   Fig.  6  shows  that  the  NVA  height  changes  necessary 
to  generate  gradients  consistent  with  the  winds  are  both 
positive  and  negative,  and  the  changes  are  similar  in 
scale  to  the  wind  maxima  themselves.   By  contrast,  the 
changes  in  the  two  cases  using  the  nonlinear  balance  equa- 
tion approach  are  on  a  much  larger  scale  as  shown  for  the 
no-flux  case  in  Fig.  7.   The  patterns  are  on  the  order  of 
the  long  waves  that  are  present,  and  the  changes  tend  to 
deepen  the  troughs  by  building  a  ridge  downwind  of  the 
trough  axis.   Unlike  the  NVA  results,  all  height  changes  are 
positive  in  the  northern  hemisphere  and  are  about  four  times 


31 


Fig.  6.   Height  changes  resulting  from  NVA  initialization 
program  (1,7,25).   Units  =  m. 


Fig.  7.   Height  changes  resulting  from  nonlinear  balance 
solution  with  no-flux  boundary  conditions. 
Units  =  m. 
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as  large  as  those  in  the  NVA.   Area  average  changes  are 
slightly  larger  in  the  no-flux  boundary  condition  than  in 
the  restoration  case,  since  the  boundary  heights  are  con- 
stant . 

Early  development  experiments  made  use  of  climatologi- 
cal  data  as  mentioned  before.   One  would  expect  these  fields 
to  be  composed  of  long  wave  features  with  none  of  the  short 
wave  systems  apparent  in  daily  analyses.   This  is  evident 
in  the  wind  and  height  changes  necessary  to  balance  the 
real  data  fields.   The  wind  changes  due  to  initialization 
of  the  climatological  data  are  one-fourth  to  one-half  of 
those  for  real  data. 

The  second  objective  of  the  diagnostic  experiments  was 
to  develop  weighting  values  for  the  variational  analysis 
scheme  which  produced  fields  with  sufficient  balancing  to 
successfully  initialize  the  prognostic  program.   Table  2 
shows  the  results  of  this  experimentation.   Since  only  a 
wind  analysis  is  available  in  the  tropics,  the  NVA  analysis 
should  be  forced  toward  the  wind  values  with  little  weight 
on  the  height  values.   The  dynamical  constraint  should  be 
large  enough  to  assure  compatibility  as  an  initial  state 
for  the  prediction  model.   One  must  pay  the  price  of  addi- 
tional computer  time  if  the  analysis  is  to  be  forced  toward 
exact  balance,  or  small  deviations  from  the  input  wind 
fields,  as  shown  in  Fig.  8.   An  arbitrary  root-mean-square 
deviation  in  the  wind  field  of  1.0  m/sec  was  set  as  a 
desirable  limit  on  changes  allowed  from  the  input  value. 
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TABLE  2 

VARIATIONAL  ANALYSIS  WEIGHTS 

a  =  winds 

3  =  geopotential 

a  =  dynamics 

Balancing  Area     Area     RMS      Area     Area     RMS 
Weights    Average  Average  Wind     Average  Average  Height 
0,  a,  a    Wind     Modulus  Change   Height   Modulus  Change 
Change   Wind     (m/sec)  Change   Height   (m) 
(m/sec)  Change           (m)      Change 
(m/sec) (m) 

1 
1,1^\25    u=1.43   1.40     2.47     -16.5    29.0     35.6 
v=  .13 

1,7,25     u=  .52    .63     1.10     -22.1    31.0     37.0 
v=  .025 

1,16,25    u=  .26    .33      .60     -25.2    34.8     42.7 
v=  .03 

1,25,25    u=  .17    .22      .41     -26.5    35.6     43.7 
v=  .02 
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Fig.  8.   NVA  wind  weighting  factor  (solid),  and  number  of 
iterations  required  for  convergence  (dashed)  for 
a  RMS  wind  difference  between  final  and  initial 
fields.   The  values  of  pressure  and  dynamical 
weighting  factors  are  held  fixed. 
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This  resulted  in  an  a  value  of  approximately  7.0  when  com- 
bined with  3  =  1.0  and  a  =  25.0  where  3  limits  height 
changes,  a  specifies  the  degree  of  horizontal  acceleration 
in  the  wind  field  and  a  limits  the  wind  changes  as  speci- 
fied in  Eqn.  9.   One  should  note  that  limiting  the  wind 
deviations  to  very  small  changes  causes  a  rapid  increase  in 
the  number  of  iterations  required  for  convergence.   Speci- 
fying the  coefficients  of  1,7,25  required  approximately  23 
seconds  of  computer  time  for  each  iteration  of  the  uncon- 
verged  points.   This  results  in  a  total  run  time  of  about 
32  minutes  which  is  extremely  long  for  an  operational  analy- 
sis.  In  contrast,  the  nonlinear  balance  approach  requires 
approximately  eight  minutes  of  computer  time. 

It  was  discovered  in  the  experiments  with  climatologi- 
cal  data  that  the  wind  fields  resulting  from  the  variational 
analysis  program  produced  excessive  two-gridlength  gravity 
waves  and  excessive  vertical  motions.   Eventually,  this 
noise  led  to  computational  instability.   The  NVA  procedure 
satisfies  the  dynamical  constraints  on  the  fields  by  limit- 
ing the  magnitude  of  horizontal  accelerations,  which  are  ob- 
tained from  the  wind  tendency  equations.   There  is  no  direct 
control  on  the  amount  of  divergence  allowed.   Therefore,  the 
third  objective  of  the  initialization  experiments  was  to 
evaluate  various  methods  of  limiting  divergence  in  the  NVA 
output  field.   Two  alternatives  existed;  make  the  input 
wind  fields  to  the  NVA  analysis  nondivergent  and  produce  a 
divergent  output;  or  make  the  output  of  the  NVA  program  non- 
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divergent  by  using  only  the  nondivergent  component.   However 
this  pos tad jus tment  of  the  wind  will  create  some  inconsist- 
ency with  the  height  field  that  was  produced  by  the  NVA 
scheme . 

A  FNWC  subroutine,  which  uses  version  III  of  the  method 
described  by  Hawkins  and  Rosenthal  (1965),  was  used  to  com- 
pute the  nondivergent  wind.   It  should  be  noted  that  this 
approach  does  not  retain  the  original  u  and  v  winds  on  the 
boundary  as  boundary  conditions.   Nevertheless,  changes  are 
minimal  on  the  boundaries. 

Table  3  shows  the  results  of  these  experiments.   The 
patterns  of  wind  changes  necessary  to  make  the  fields  non- 
divergent  show  patterns  similar  to  those  of  the  nonlinear 
balance  equation,  however,  on  a  smaller  scale,  as  shown  in 
Fig.  9.   The  individual  wind  maxima  are  altered  similar  to 
the  scale  of  the  NVA  changes.   One  can  see  that  the  NVA  pro- 
gram changes  the  winds  approximately  the  same  amount  whether 
the  input  is  divergent  or  nondivergent;  however,  the  height 
changes  necessary  with  the  nondivergent  input  are  larger. 
The  modulus  and  root-mean-square  wind  changes  necessary  to 
make  the  fields  nondivergent  are  about  two  and  one-half 
times  larger  than  changes  in  the  NVA  analysis.   Neverthe- 
less, the  total  change  involved  in  both  methods  are  very 
nearly  equal.   It  should  be  pointed  out  that  a  property  of 
a  nondivergent  field  is  that  a  cyclic  row  average  of  the  v 
component  wind  around  the  global  band  equals  zero.   This 
implies  no  net  mass  transport  in  the  meridional  direction. 
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TABLE  3 


VARIATIONAL  ANALYSIS  METHODS 

Method     Area     Area     RMS      Area     Area     RMS 

Average  Average  Wind     Average  Average  Height 
Wind     Modulus  Change   Height   Modulus  Change 
Change   Wind     (m/sec)  Change   Height   (m) 
(m/sec)  Change  (m)      Change 

(m/sec) (m) 

Method  1 

a. 

Divergent  u=.25    1.75     2.54     None     None     None 

Input      v= . 2  7 

Made 

Nondivergent 

b. 

Nondivergent 

Input  To   u=.50     .59     1.00     -31.4    40.2     50.8 

NVA        v=-.02 

(1,7,25) 

Total      u=.75    1.80     2.54     -31.4    40.2     50.8 
Change     v= . 25 

************************** 

Method  2 

a . 

Divergent 

Input  To   u=.52     .60     1.10     -22.0    31.0     37.0 

NVA        v=.025 

(1,7,25) 

b. 

NVA 

(1,7,25)   u=-.02   1.67     2.44     None     None     None 

Output     v=.25 

Made 

Nondivergent 

Total      u=.50    1.81     2.61     -22.0    31.0     37.0 
Change     v=.27 
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Fig.  9.   V  component  wind  speed  changes  resulting  from 
nondivergent  subroutine.   Positive  (+)  value 
indicates  increase  in  southerly  wind  speed. 
Units  =  m/sec . 
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VI.   PREDICTION  RESULTS 

A.   BOUNDARY  CONDITIONS 

Boundary  conditions  on  a  global  band  pose  problems  dif- 
ferent than  those  in  a  hemispheric  model.   The  global  band 
boundaries  are  in  the  mid-latitudes  where  strong  meridional 
flow  causes  strong  cross-boundary  fluxes.   A  hemispheric 
model  avoids  these  difficulties,  since  boundaries  are  loca- 
ted in  the  tropics  where  pressure  gradients  are  relatively 
small  and  large  scale  waves  are  not  commonly  present.   The 
major  problem  in  the  prediction  models  caused  by  the  large 
cross-boundary  flow  is  the  generation  of  gravity  waves  at 
and  near  the  boundaries.   These  waves  contaminate  the  pre- 
diction of  the  meteorological  waves  and  may  lead  to  computa- 
tional instability.   Prediction  on  the  full  grid  with  no- 
flux  boundary  conditions  produce  only  small  amplitude  gravi- 
ty waves.   Restoration  boundary  conditions  with  compatible 
nonlinear  balance  equation  input  to  the  prediction  model 
produce  strong  two-gr idlength  gravity  waves,  which  require 
filtering  to  obtain  a  stable  24-hour  forecast.   This  filter- 
ing is  accomplished  using  a  scheme  devised  by  Shapiro  (1971) 
in  which  the  variable's  value  at  J, I  is  altered  by  1/8 

<sj+i,i+  sj-i,i+  sj,i+i+  sj,i-i-  A  sj,i>  where  s  ls  a 

general  variable.   In  order  to  obtain  a  24-hour  forecast,  it 
was  necessary  to  apply  the  filter  once  every  four  hours. 
Figs.  10  and  11  show  the  gravity  wave  effect  in  the  RMS  <f> 
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Fig.  10.  <J>RMS  vs  •  Grid  Row  for  full  grid  restoration 
prediction.   Units  =  m . 
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Fig.  11.  ooRms  vs«  Grid  Row  for  full  grid  restoration 
prediction.   Units  =  mb/sec. 
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and  RMS  w  values  respectively.   The  2-d  nature  of  the  gravi- 
ty waves  is  evident  in  both  fields.   The  cross-boundary 
flow  and  resulting  gravity  waves  produce  much  larger  RMS  <f> 
and  w  values  near  the  northern  wall  which  increases  in 
amplitude  with  time.   NVA  input  to  the  full  prediction  model 
resulted  in  the  most  severe  2-d  waves,  as  might  be  expected 
since  the  dynamic  constraint  is  not  as  strong  as  in  the 
nonlinear  balance  approach.   It  was  necessary  to  apply  the 
filter  twice  at  four-hour  intervals  to  produce  a  24-hour 
forecast.   The  staggered  grid  predictions  did  not  require 
any  filtering  to  produce  a  stable  24-hour  forecast.   RMS  <J) 
and  RM S  oo  for  the  staggered  grid  restoration  prediction, 
shown  in  Figs.  12  and  13,  are  much  better  behaved  than  in 
the  full  grid  prediction.   At  t  =  0,  the  to  curve  in  Fig.  13 
arises  from  extracting  alternate  u  and  v  component  winds 
from  the  balanced  solution  obtained  on  the  full  grid.   There' 
fore,  the  staggered  field  contains  divergence  while  the  full 
grid  field  is  nondivergent  and  contains  no  vertical  motion 
at  t=0 .   The  vertical  motions  in  the  staggered  prediction 
maintain  the  same  size  03  as  time  increases.   The  amplitude 
of  the  two-gridlength  waves  is  not  as  great  and  values  at 
the  northern  boundary  do  not  increase  in  amplitude  as  rapid- 
ly with  time  as  in  the  full  grid  prediction  models.   This 
is  to  be  expected  since  the  finite  differencing  involves  cal' 
culations  over  two  gridlengths  of  the  full  grid  and  results 
in  the  use  of  points  which  have  a  similar  phase  angle  with 
respect  to  the  2-d  waves. 
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Fig.  12.  $rmS  vs.  Grid  Row  for  staggered  grid  restoration 
prediction.   Units  =  m . 
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Fig.  13.  WrMc  vs.  Grid  Row  for  staggered  grid  restoration 


prediction.   Units  =  mb/sec. 
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B.   PREDICTION  VERIFICATION  METHODS 

The  output  fields  were  smoothed  with  a  filter  which 

altered  the  parameter  value  at  J, I  by  1/8  (S   .  T+S   , 

+  S T  t.t+S       -4  S  „  T).   Verification  of  the  prognosis 
J.I+1  Jtl_±  J, I 

was  based  on  these  filtered  fields.   Also,  to  avoid  un- 
realistic fields  produced  by  the  restoration  process  near 
the  boundaries,  and  fictitious  boundary  specification  in 
the  objective  analysis,  verification  statistics  were  only 
computed  in  the  interior  region  where  a  true  dynamical 
forecast  was  computed.   This  excluded  the  outer  six  rows 
in  the  full  grid  prediction  and  the  outer  four  rows  in  the 
staggered  grid.   Verification  criteria  for  predictive 
models  are  numerous  and  by  no  means  obvious.   The  area 
average  modulus  wind  component  differences  and  the  RMS  dif- 
ferences were  chosen  since  accuracy  in  the  wind  analysis  is 
more  desirable  in  the  tropics  than  are  good  height  values. 
A  single  numerical  value  which  was  used  to  combine  the 
various  measures  of  forecast  ability  was  defined  as  in  the 
following  formula: 

Forecast  Value  =  ( I u I „     ■   -  I u I _   ,  _)  + 

■  ' Persist   '  ' Verxf 

(lvlpersist-lvlVerif)  +  <U         "  U      } 

verii       RMg         RMS 

Persist     Verif 

+  (v        -  v      )  (10) 

RMS         RMS 
Persist    Verif 

where  |u|  =  area  average  modulus  error 

URMS  ~  root  mean  square  error 

Persistence  forecast  error  is  defined  as  the  comparable 
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balanced  field  at  t+24  hours  minus  the  comparable  balanced 

o 

field  at  t  .   Lavoie  and  Wiederanders  (1968)  have  pointed 
out  that  in  the  tropics  a  persistence  forecast  plus  a  varia- 
ble amount  of  climatology  is  still  difficult  to  beat  with 
any  objective  prediction  technique. 

Typical  magnitude  persistence  forecast  error  values  are 
shown  in  Table  4. 

Verifying  forecast  error  is  equal  to  the  predicted 

values  at  t  +24  hours  minus  the  verifying  balanced  field  at 
o 

t+24  hours.  Therefore,  if  the  forecast  value  is  positive, 
the  verification  is  better  than  persistence  with  increasing 
magnitudes  denoting  better  verification  skill. 

A  negative  value  implies  less  skill  than  a  persistence 
forecast.   Less  emphasis  was  placed  on  height  verification 
since  these  fields  are  derived  rather  than  analysed.   How- 
ever, one  measure  of  prediction  accuracy  is  the  ratio  of 
RMS  persistence  height  change  to  the  RMS  verification 
height  error.   Numbers  less  than  one  imply  a  persistence 
forecast  is  more  accurate.   Values  of  this  ratio  greater 
than  one  imply  forecast  skill  compared  to  persistence. 
These  values  are  shown  in  Tables  5  and  6  in  parentheses. 

The  number  of  passes  through  the  predicted  fields  with 
the  smoother  had  a  definite  effect  on  the  verification 
statistics.   A  typical  example  with  the  forecasts  with  NVA 
input  is  shown  in  Table  4.   The  modulus  and  RMS  verification 
values  for  the  v  component  are  better  than  persistence, 
regardless  of  the  number  of  filter  passes.   Successive 
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TABLE  4 

VERIFICATION  STATISTICS  AS  FUNCTION  OF  FILTER  PASSES 
Forecast:  Non-divergent  NVA  input  (1,7,25) 

Number  Area  Area     RMS  RMS        Area      RMS 

of     Average  Average  u  v  Average   Height 

Passes  Modulus  Modulus  Component  Component  Modulus   Differ- 
u  v        Wind  Wind       Height    ence 
Compon-  Compon-  Dif f erenceDif f erenceDif f er-   (m) 
ent  ent      (m/sec)  (m/sec)    ence 
Wind  Wind  (m) 
Differ-  Differ- 
ence ence 
(m/sec)  (m/sec) 


Persis 

tence 

4.19 

5.45 

7.01 

10.25 

31.35 

63.40 

None 

4.82 

4.35 

7.05 

7.25 

51.05 

71.66 

1 

4.46 

4.15 

6.48 

6.94 

49.29 

68.89 

2 

4.32 

4.07 

6.24 

6.78 

48.57 

67.73 

3 

4.27 

4.04 

6.13 

6.70 

48.18 

67.09 

4 

4.27 

4.02 

6.10 

6.65 

47.97 

66.  70 

5 

4.30 

4.03 

6.12 

6.63 

47.92 

66.60 

6 

4.34 

4.04 

6.16 

6.64 

48.02 

66.66 

7 

4.40 

4.04 

6.23 

6.64 

48.15 

66.75 
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filtering  up  to  four  or  five  times  improves  the  verifica- 
tion values.   However,  beyond  five  filterings,  the  verifi- 
cation values  begin  to  decline.   The  u  wind  component  veri- 
fication is  different,  however.   Area  average  modulus  u 
component  differences  remainhigher  than  persistence  values. 
RMS  u  component  differences  with  one  pass  of  the  filter  are 
smaller  than  for  a  persistence  forecast  and  continue  to 
improve  up  to  four  or  five  passes.   This  can  be  explained  by 
noting  that  wind  gradients  north  and  south  of  u  component 
wind  maximum  are  much  larger  than  the  gradients  around  v 
component  wind  maximum  (see  Figs.  14  and  16).   Therefore, 
a  slight  displacement  of  the  predicted  maximum  north  or 
south  of  the  actual  position  leads  to  larger  errors  than  are 
produced  in  the  v  component  by  east-west  phase  errors. 
Even  though  filtering  reduces  the  actual  jet  core  speed, 
the  diffusion  of  the  maximum  due  to  filtering  results  in  a 
better  verification  value. 

Figs.  14  through  19  show  the  fields  used  for  NVA  prog- 
nosis verification,  and  the  prognosis  error  fields  for  the 
v  component,  u  component  and  height  fields  respectively. 
All  fields,  especially  the  v  component  and  height  field 
show  that  horizontal  advection  of  most  of  the  synoptic 
features  was  too  rapid.   Predicted  positions  and  intensity 
of  the  major  synoptic  features  are  shown  by  the  stars  and 
values  in  parentheses  on  the  verification  charts.   One  can 
see  that  the  analysed  error  values  are  predominantly  a 
function  of  position,  and  not  intensity  or  addition/deletion 
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Fig.  14.  300  mb  NVA  balanced  v  component  wind  verification 
analysis  from  60E  to  150W  for  09/1200Z  January 
1972.   Star  denotes  prognosis  wind  maxima  location 
with  intensity  in  parenthesis.   Units  =  m/sec. 


Fig.  15.  Nondivergent  input  NVA  prediction  verification 
error  of  v  component  wind.   Positive  (+)  value 
indicates  too  strong  southerly  wind.   Units  = 
m/sec  . 
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Fig.  16.  300  mb  NVA  balanced  u  component  wind  verification 
analysis  from  60E  to  150W  for  09/1200Z  January 
1972.   Star  denotes  prognosis  wind  maxima  location 
with  intensity  in  parenthesis.   Units  =  m/sec. 
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Fig.  17.  Nondivergent  input  NVA  prediction  verification 
error  of  u  component  wind.   Positive  (+)  value 
indicates  too  strong  westerly  wind.   Units  = 
m/sec . 
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Fig.  18.  300  mb  NVA  balanced  height  verification  field  from 
60E  to  150W  for  09/1200Z  January  1972.   Trough 
symbol  denotes  prognosis  trough  position.   Units  = 
decame  ters . 
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Fig.  19.  Nondivergent  input  NVA  prediction  verification 
error  of  height  field.   Units  =  m. 
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of  systems.   A  typical  example  is  the  short  wave  trough 
southwest  of  Korea.   The  prediction  model  overpredicts  the 
trough  movement  (See  Fig.  18)  with  a  resulting  positive 
height  error  pattern  west  of  the  predicted  trough  position 
and  a  negative  center  east  of  the  trough  (see  Fig.  19). 
The  model  moves  both  of  the  v  component  maxima  associated 
with  the  trough  too  far  south  as  shown  in  Fig.  14.   This 
position  error  contributes  about  one-half  of  the  value  of 
the  verification  error  analysis  shown  in  Fig.  15.   The  re- 
mainder of  the  error  is  due  to  an  underestimate  of  the 
magnitude  of  the  v  component  maxima,  which  seems  to  be  true 
of  the  majority  of  the  centers.   The  u  component  wind 
minima  associated  with  the  trough  is  forecast  to  move  more 
rapidly  eastward  than  the  verifying  position  shown  in  Fig. 
16.   This  is  depicted  in  Fig.  17  as  a  negative  error  center 
east  of  Korea  and  a  positive  center  west  of  Korea.   Unlike 
the  v  component  winds,  the  prediction  appears  to  retain  u 
component  magnitudes  comparable  with  the  verification. 

C.   RESTORATION  TECHNIQUES 

Results  of  the  experiments  with  climatological  input 
data  had  suggested  that  restoration  of  both  wind  fields  and 
the  height  field  at  each  time  step  also  tended  to  generate 
excessive  2-d  gravity  waves  originating  from  the  boundaries 
This  can  be  explained  as  follows:   As  the  wind  and  height 
fields  in  the  region  of  outflow  try  to  adjust  to  the  fixed 
boundary  values,  the  restoration  coefficient  allows  only 
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partial  adjustments.   This  produces  a  continual  imbalance, 
which  is  manifested  by  gravity  waves  propagating  from  this 
region.   The  generation  of  these  unwanted  gravity  waves 
can  be  reduced  by  either  restoring  the  wind  values  and  not 
the  heights,  or  vice-versa.   Slightly  smaller  amplitude 
gravity  waves  were  produced,  and  better  verification  values 
were  obtained,  when  u  and  v  component  winds  were  restored 
with  no  restoration  of  the  height  field.   This  result  is 
shown  in  Table  5A  for  both  the  full  and  staggered  grid  re- 
sults where  Forecast  value  is  defined  as  in  Eq.  10  and  the 
value  in  parenthesis  is  the  ratio  of  RMS  height  error  to  a 
RMS  height  error  from  a  persistence  forecast. 

D.   DIVERGENCE  VALUES  OF  INPUT  NVA  FIELDS 

As  mentioned  earlier,  a  limitation  of  the  divergence  in 
the  NVA  input  to  the  prediction  model  appeared  to  have  de- 
sirable effects  on  the  elimination  of  unwanted  two-grid- 
length  gravity  waves.   This  might  be  expected  to  have  an 
effect  on  the  verification  values  of  the  prediction,  and  is 
shown  in  Table  5B.   The  best  verification  was  obtained  using 
an  NVA  output  made  nondivergent  input  to  the  prediction 
(Method  2  in  Table  3).   Slightly  worse  results  were  obtained 
from  NVA  input  data  obtained  from  a  nondivergent  input  into 
the  NVA  program  (Method  1  in  Table  3).   Considerably  less 
skill  was  achieved  using  initial  data  from  the  NVA  program 
with  divergent  input  (Method  2a  in  Table  3).   In  fact,  the 
full  grid  prediction  model  blew  up  before  a  24-hour  forecast 
could  be  obtained.   One  would  expect  a  wind  field  which  was 
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TABLE  5 


PROGNOSIS  VERIFICATION  AS  A  FUNCTION  OF 
RESTORATION,  NVA  DIVERGENCE,  NVA  WEIGHTING 


Input 


Restore 


Forecas  t 
Value 


A.  RESTORATION 

Full  Grid 

Restoration 
Restoration 

Staggered  Grid 

Restoration 
Restoration 
Restoration 

NVA(1 , 7 , 25 )made  nondivergent 
NVA(1 , 7 , 25 ) made  nondivergent 
NVA(1 , 7 , 25 )raade  nondivergent 

B.  NVA  DIVERGENCE 

Full  Grid 

NVA(1 , 7 , 25 ) wi th  nondivergent  input 
NVA(1 , 7 , 25) made  nondivergent 

Staggered  Grid 

NVA( 1 , 7 , 25) wi th  nondivergent  input 
NVA(1 , 7 , 25) made  nondivergent 

C.  NVA  WEIGHTING 

Full  Grid 
NVA(1, 1-1/2, 25) 

NVA(1,7,25) 

NVA(1,25,25) 

Staggered  Grid 
NVA(1, 1-1/2, 25) 
NVA(1,7,25 
NVA(1,25,25) 


u  &  v 


u  &  V 
U  &  V 


U  &  V 
U  &  V 


U  &  V 
U  &  V 
U  &  V 


U  &  V 
U   &  V 

u  &  v 


2.08  + 
2.13  + 


37 
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+ 
+ 


(.77) 
(.68) 


All 
u     &     V 

2.82  + 
2.90  + 
3.77    + 

(.61) 
(.65) 
(.58) 

4> 

All 
u     &    V 

5.69  + 
5.80  + 
5  .86    + 

(.91) 
(.93) 
(.95) 

(.69) 
(.84) 


73  +  (.96) 
86  +  (.95) 


Computational 
Ins  tability 
Computat  ional 
Ins  tability 
3.30  +  (.49) 


3.34  +  (.78) 
3.57  +  (.75) 
3.91  +  (.72) 
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nondivergent  to  produce  fewer  gravity  waves  than  a  divergent 
wind  field,  and  the  resulting  gravity  waves  appear  to  be  one 
factor  responsible  for  degrading  the  verification  results  in 
the  divergent  NVA  input  wind  fields. 

E.  EFFECT  OF  NVA  WEIGHTING 

A  variation  of  the  weighting  factors  in  the  NVA  analysis 
could  be  expected  to  vary  the  prognosis  verification.   A 
closer  specification  of  the  wind  changes  allowed  in  the  NVA 
analysis  (larger  a)  should  produce  a  better  forecast  wind 
field,  as  described  in  the  diagnostic  results  section  on 
NVA  weighting.   Table  5c  shows  this  result  for  the  stag- 
gered grid  prediction.   In  the  full  grid,  specifying  5=25 
produced  a  computationally  stable  prediction  for  24-hours, 
but  using  B  =  7  or  1.5  did  not  adequately  limit  the  changes 
in  the  wind  components  and  the  predictions  became  unstable 
before  24-hours. 

F.  STAGGERED  VS.  NON-STAGGERED  GRID 

Tables  5  and  6  point  out  that  the  staggered  grid  predic- 
tions in  all  cases  produce  more  accurate  forecasts  than  a 
full  grid  with  comparable  input.   It  is  difficult  to  deter- 
mine if  this  is  a  result  of  the  removal  of  the  computational 
mode  from  the  walls  by  use  of  the  staggered  grid  or  an  even 
more  severe  overprediction  of  the  movement  of  systems  in  the 
full  grid  compared  to  the  staggered  grid.   The  excessive 
amplitude  of  the  2-d  gravity  waves  evident  in  the  full  grid 
predictions,  as  compared  to  the  staggered  grid,  appears  to 
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be  the  most  significant  influence. 

G.   COMPARISON  OF  INITIALIZATION  METHODS 

On  both  the  staggered  and  the  nonstaggered  grids,  the 
NVA  fields  made  nondivergent  input  to  the  prediction  had 
the  best  verification  value.   Considerably  poorer  were  the 
balance  restoration  and  balance  no-flux,  with  the  forecasts 
with  restoration  being  slightly  better  (see  Table  6) . 
Better  verification  of  the  NVA  input  might  be  expected 
since  the  changes  necessary  in  the  diagnostic  phase  of 
prediction  are  smaller  for  winds  and  much  smaller  for 
height  than  in  the  nonlinear  balance  method.   This  results 
in  individual  synoptic  features  retaining  their  identity 
and  as  a  result  are  forecast  more  accurately. 
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Input 


Full  Grid 


TABLE  6 


PROGNOSIS  VERIFICATION  AS  A  FUNCTION  OF 
INITIALIZATION  TYPE 


Restore    Forecast 
Value 


No-Flux 

Restoration 

NVA(1 , 7 , 25 )made  nondivergent 


.20  +  (.57) 
u  &  v  2.13  +  (.68) 
u  &  v      4 .52  +  ( . 84) 


Staggered  Grid 


No-Flux 

Restoration 

NVA(1 , 7 , 25) made  nondivergent 


3.55  +  (.55) 
u  &  v      3.77   +  (.58) 
u  &  v      5 . 86  +  (.95) 
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VII.   CONCLUSIONS  AND  RECOMMENDATIONS 

The  staggered  grid  forecasts  produced  better  verifica- 
tion statistics  than  the  full  grid.   This  is  desirable  since 
the  barotropic  full  grid  prediction  requires  approximately 
33  minutes  of  computer  time  for  a  24-hour  forecast,  while 
the  staggered  grid  requires  only  about  10  minutes.   This  is 
a  considerable  saving  of  time  when  considering  operational 
requirements . 

Restoration  near  the  boundaries  of  the  u  and  v  component 
winds  and  no  restoration  of  the  height  field  produced  better 
verification  values  than  restoring  all  fields  or  restoring 
only  the  height  field.   This  results  in  a  slight  additional 
reduction  in  computer  time. 

A  nondivergent  NVA  input  to  the  prediction  model  pro- 
duced better  verification  statistics  than  either  the  balance 
restoration  or  balance  no-flux  technique.   Nondivergent  NVA 
input  also  produced  better  verification  statistics  than  NVA 
fields  which  had  no  constraint  on  divergence.   The  smallest 
amount  of  adjusting  of  wind  speed  from  the  input  fields, 
offset  by  rapidly  increasing  NVA  convergence  time,  produced 
better  verification  values.   To  obtain  these  NVA  fields  re- 
quired running  two  programs  with  a  combined  run  time  of 
about  thirty-five  minutes.   Again,  this  is  excessive  for  an 
operational  product.   Three  changes  are  possible  to  reduce 
the  computation  time.   First,  the  NVA  program  could  be 
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altered  to  place  an  additional  constraint  on  the  magnitude 
of  the  divergence  of  the  resulting  fields.   This  would 
eliminate  the  need  for  running  the  nondivergence  program. 
Also,  further  experiments  might  show  that  a  specified 
amount  of  divergence  results  in  a  better  initial  field  than 
no  divergence.   Also,  time  might  be  reduced  by  altering  the 
values  of  the  NVA  coefficients  to  require  less  adherence  to 
the  input  fields,  and  possibly  larger  tolerances  on  the 
horizontal  acceleration.   The  most  significant  time  saving 
could  be  realized  by  converting  the  NVA  program  to  operate 
on  a  staggered  grid.   This  would  allow  for  a  solution  con- 
sidering only  one-quarter  of  the  present  points.   Since 
iterative  processes  are  inherently  more  efficient  on  a 
smaller  grid,  the  time  savings  should  be  well  in  excess  of 
a  factor  of  four. 

Boundary  conditions  on  the  Mercator  belt  still  remain 
a  problem.   Research  by  Lieutenant  Commander  E.  Harrison 
(personal  communication)  on  appropriate  boundary  conditions 
for  numerical  models  has  suggested  a  different  method  of 
handling  the  prediction  in  regions  of  outflow. 

It  is  to  be  expected  that  a  barotropic  model,  which  can 
only  advect  synoptic  features  with  the  wind  at  one  level, 
should  produce  errors  in  displacement;  specifically  a  300  mb 
level  should  overforecast  movement.   In  the  real  atmosphere 
the  vertical  differences  in  both  temperature  and  momentum 
advection  can  make  the  barotropic  approximation  a  poor  one. 
This  would  suggest  that  a  baroclinic  model  would  give 
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significantly  better  movement  forecasts.   The  staggered  grid 
now  makes  a  three-  or  four-level  model  operationally  feasi- 
ble.  Initial  experiments,  on  a  limited  grid  used  by  Harri- 
son (1969)  and  Steinbruck  (1971),  have  shown  that  at  least 
a  three-level  model  is  desirable.   At  present,  the  time  step 
in  the  staggered  prediction  models  is  five  minutes.   This 
could  be  doubled  if  the  advection  of  the  height  field  in 
Eq.  5  was  finite-differenced  over  twice  the  present  space 
increment  in  the  full  grid,  as  is  the  case  for  momentum  ad- 
vection, with  a  staggered  grid  in  the  time  sense.   The  use 
of  semi-implicit  time  differencing  in  the  prediction  models 
would  also  be  a  means  for  increasing  the  time  step.   One 
possible  approach  is  a  time-averaging  form  suggested  by 
Brown  (1971). 

By  using  at  least  three  layers,  it  is  possible  to 
parameterize  the  effect  of  cumulus  convection,  which  is  so 
important  in  the  tropics.   A  baroclinic  model  would  also 
make  experiments  with  tropical  cyclones  movement  and 
development  feasible. 
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